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Abstract 

The  least-squares  analysis  of  Two-Way  Satellite  Time  and  Frequency  (TWSTFT)  data,  with 
the  aim  of  determining  optimal  estimates  of  phase  and  frequency  offsets,  is  considered.  An 
overview  of  the  application  of  Gauss-Markov  estimation  to  the  analysis  of  a  uniformly  spaced 
time  series  is  presented.  Two  aspects  of  TWSTFT  data  analysis  are  examined  in  depth.  Firstly, 
measures  based  on  second  difference  statistics  for  characterizing  the  consistently  but  unevenly 
spaced  TWSTFT  measurements  are  introduced.  Secondly,  an  approach  to  the  estimation  of  the 
occasional  and  unknown  delay  steps  in  TWSTFT  data,  for  example  due  to  hardware 
replacement,  is  presented. 


1  INTRODUCTION 

Two-Way  Satellite  Time  and  Frequency  Transfer  (TWSTFT)  [1]  has  been  used  operationally  by  primary 
timing  laboratories  for  several  years.  TWSTFT  links  between  laboratories  have  proved  to  be  both  stable 
and  reliable,  and  consequently  TWSTFT  measurements  are  used  for  several  of  the  main  links  in  the 
computation  of  International  Atomic  Time  (TAI).  This  paper  describes  new  analysis  techniques  being 
developed  for  the  processing  of  TWSTFT  measurements.  The  overall  aim  of  the  study  is  to  be  able  to 
determine,  at  any  measurement  epoch,  estimates  of  the  phase  and  normalized  frequency  offsets,  together 
with  their  uncertainties,  between  two  clocks  or  time  scales  being  compared  using  TWSTFT. 

In  Section  2,  we  describe  an  approach,  based  on  Gauss-Markov  estimation,  to  solving  this  problem  for  the 
case  of  measurements  of  the  phase  difference  between  two  clocks  or  time  scales  made  at  uniformly  spaced 
times  and  assuming  knowledge  of  the  random  noise  processes  underlying  the  measurements.  For  time 
series  corresponding  to  uniformly  spaced  times,  second  difference  statistics,  such  as  the  Allan  and 
modified  Allan  variances,  provide  a  means  to  obtain  information  about  the  noise  processes. 

However,  the  application  of  this  approach  to  TWSTFT  is  made  difficult  for  two  reasons.  Firstly,  for  many 
operational  links  based  on  TWSTFT,  measurements  are  available  only  on  Mondays,  Wednesdays,  and 
Fridays.  The  measurements  constitute  a  time  series  that  is  unevenly  spaced  (with  a  2-,  2-  and  3-day 
spacing),  but  consistent,  in  that  this  spacing  pattern  is  repeated.  Conventional  second  difference  statistics, 
such  as  the  Allan  and  modified  Allan  variances,  are  not  applicable  to  such  time  series.  A  new  second- 
difference  statistic,  designed  for  consistently  but  unevenly  spaced  time  series,  is  described  in  Section  3. 
Secondly,  it  is  well  known  that  TWSTFT  data  can  be  subject  to  occasional  phase  offsets  or  steps  resulting 
from  replacement  of  the  hardware  underlying  the  collection  of  the  data.  In  Section  4,  a  method  is 
described  for  analyzing  TWSTFT  measurements  to  provide  estimates  of  these  steps,  together  with  their 
associated  uncertainties.  A  summary  is  given  in  Section  5. 
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Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 
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2  ANALYSIS  OF  A  UNIFORMLY  SPACED  TIME  SERIES 


Let  x(t)  denote  the  phase  difference,  as  a  function  of  time  t,  between  two  clocks  and  {xt:  i  =  1,  m]  a 
time  series  of  measurements  of  x  made  at  uniformly-spaced  times  {t,:  i  =  1,  m}.  Suppose  x(t)  is 

modeled  in  terms  of  a  phase  offset  p0  and  a  normalized  frequency  offset  /0,  i.e., 

x(t)  =  p0+f0t,  (1) 

and  the  measurements  as 

=  Po+/o*,+<v  (2) 

where  e,  denotes  a  random  error.  Using  matrix  notation, 

x  =  Xp+e,  (3) 

where 


M 

f  l 

O 

X  = 

*2 

,  x  = 

l 

h 

,  p  = 

W 

,  e  = 

*2 

\XmJ 

) 

\foj 

KemJ 

We  suppose  that  the  elements  of  e  have  zero  expectation  (E(e)  =  0)  and  a  known  covariance  matrix  (V(e)  = 
V)  [2].  Typically,  the  elements  e,  are  a  linear  combination  of  (independent)  samples  of  known  noise  types. 
For  simplicity  of  presentation,  we  restrict  the  analysis  to  the  well-known  noise  types  of  white  phase 
modulation  (WPM),  white  frequency  modulation  (WFM),  and  random-walk  frequency  modulation 
(RWFM).  The  extension  of  the  work  to  other  noise  types,  such  as  flicker  phase  modulation  (FPM)  and 
flicker  frequency  modulation  (FFM),  given  a  characterization  of  these  noise  types,  is  straightforward. 

We  express  V  in  terms  of  known  covariance  matrices  that  relate  to  each  noise  type  and  parameters 
describing  the  ‘‘magnitude”  of  each  noise  type  present.  Each  noise  type  is  characterized  by  a 
transformation  of  a  WPM  process.  The  covariance  matrices  V  for  standardized  WPM,  WFM,  and  RWFM, 
for  which  the  transformed  WPM  process  has  unit  standard  deviation,  are  given  by,  respectively, 

^WPM=^’  ^WFM  =TT  ,  VRWFM  ={r  )  , 


where  /  is  the  identity  matrix,  and  T  defines  the  “summation  operator” 

0  ...  0^ 

a,  1  1  **.  : 

T=  :  •  •  0 

,1  ...  1  1, 

[3,  Appendix  A].  The  covariance  matrix  V(e)  for  noise  arising  as  a  linear  combination  of  (independent) 
WPM,  WFM,  and  RWFM  is  then  expressed  as 
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VXe)  -  CTWPM^WPM  +  aWFM  ^WFM  +  ^RWFM  ^RWFM  ' 
where  ct2wpm,  ct2wfm,  and  ct2rwfm  denote  the  variances  of  the  WPM  underlying  each  noise  type. 


(4) 


Given  the  model  (2)  for  the  measurements,  and  assuming  the  covariance  matrix  V(e)  is  also  given,  for 
example  as  defined  by  (4),  a  leas-squares  analysis  or  Gauss-Markov  regression  [2]  may  be  applied  to 
obtain  estimates  of  the  phase  offset  p0  and  the  normalized  frequency  offset  /0,  together  with  their  associated 
uncertainties.  The  estimates  are  optimal  in  the  sense  of  being  the  estimates  of  minimum  variance  from  the 
class  of  linear  unbiased  estimates  of  po  and  f0  [2]. 

The  results  of  this  analysis  for  the  examples  of  (pure)  WPM,  WFM,  and  RWFM  are  given  in  Table  1.  For 
ease  of  presentation  of  results,  we  suppose  tt  =  ix0  with  i  =  -m, ...,  +m  for  WPM  and  i  =  0, ...,  m  for  WFM 
and  RWFM.  Other  cases  can  also  straightforwardly  be  handled.  The  results  show  that  for  WFM  and 
RWFM,  the  estimate  po  of  the  phase  offset  at  t  =  0  is  the  measurement  x0,  whereas  for  WPM  and  data 
distributed  symmetrically  about  t  =  0,  it  is  the  mean  of  the  measured  phase  offset  values.  Furthermore,  for 
WFM  and  RWFM,  the  estimate  /0  of  the  normalized  frequency  offset  is  the  slope  of  the  chord  joining  the 
first  and  last  measured  values  (for  WFM)  and  the  first  and  second  values  (for  RWFM). 


The  law  of  propagation  of  uncertainty  [4]  is  applied  to  the  model  defined  by  the  Gauss-Markov  regression 
and  the  statistical  model  (4)  for  the  measurements  to  provide  the  standard  uncertainties  of  the  estimates  po 
and  /o  and  their  covariance.  (For  the  cases  considered  in  Table  1,  the  estimates  p0  and  /0  are  uncorrelated 
for  WPM  and  WFM,  but  they  are  correlated  for  RWFM.)  The  estimates  of  po  and  /o,  together  with  the 
model  (1),  can  be  used  to  predict  the  phase  offset  jc  at  any  given  time  t.  A  further  application  of  the  law  of 
propagation  of  uncertainty  permits  the  standard  uncertainty  of  any  quantity  derived  from  p0  and  /o,  such  as 
a  predicted  phase  offset,  to  be  evaluated. 
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Table  1:  Results  of  Gauss-Markov  estimation  for  the  examples  of  (pure)  WPM,  WFM 
and  RWFM. 


3  SECOND  DIFFERENCE  STATISTICS 

For  several  years  TWSTFT  measurements  made  using  Intelsat  satellites  between  European  primary  timing 
laboratories,  and  between  European  and  North  American  laboratories,  have  been  available  on  Mondays, 
Wednesdays,  and  Fridays.  The  measurements  constitute  a  consistently  but  unevenly  spaced  time  series. 
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Time  transfer  measurements  are  usually  characterized  through  the  use  of  the  well-known  second  difference 
statistics  AVAR,  MVAR  and  TVAR  [5].  However,  these  statistics  do  not  lend  themselves  directly  to  the 
characterization  of  unevenly  spaced  measurements.  The  characterization  of  unevenly  spaced  TWSTFT 
data  has  been  considered  previously,  using  techniques  based  on  interpolating  missing  data  followed  by  the 
application  of  second  difference  statistics  to  the  resulting  “reconstructed”  time  series  [3,  6].  The  approach 
considered  here  is  to  develop  a  second-difference  statistic  that  may  be  used  directly  on  unevenly  spaced 
TWSTFT  data,  thus  avoiding  any  dependence  on  the  form  of  interpolation  used. 


3.1  Allan  Variance  for  a  Uniformly  Spaced  Time  Series 


In  order  to  motivate  the  construction  of  a  second-difference  statistic  for  a  nonuniformly  spaced  time  series, 
we  begin  by  reviewing  the  Allan  variance  for  a  uniformly  spaced  time  series  and,  in  particular,  its 
definition  in  terms  of  a  second-difference  operator. 


Consider  the  time  series  x  =  (jq,  ...,  xm)T  with  x,  denoting  a  measurement  of  x(t)  at  t,  =  rt0,  i  =  1,  ...,  m. 
Here,  x  =  x0  denotes  the  spacing  (in  time)  between  the  measurements  x  of  phase  difference.  To  determine 
the  “single-spaced”  (x  =  x0)  Allan  variance  statistic,  we  first  form  the  time  series  of  second  differences  y  = 
(yi, ....  ym-if  of  x,  viz.. 


In  matrix  form, 
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y  =  AtX, 


where  A!  is  a  band  diagonal  matrix  with  band  defined  by 

=(l  -2  l) 


The  time  series  y  is  stationary  for  the  five  well  known  noise  types  listed  in  Section  2  as  well  as  linear 
frequency  drift  (the  latter  modeled  by  behavior  of  the  form  fDt2),  and  is  independent  of  the  phase  and 
normalized  frequency  offsets  p0  and  f0  (but  not  fD)  in  a  model  for  the  original  time  series  x.  The  Allan 
variance  AVAR(x0)  is  then  calculated  from  y  using 


AVAR(x0)  =  -^E\y2 } 

2x0 

where  E\y2]  denotes  the  sample  expectation  (or  arithmetic  mean)  of  the  elements  y2  derived  from  y. 


The  Allan  variance  for  a  longer  averaging  time  x  =  nx0  may  similarly  be  constructed  using  a  second- 
difference  operator  a„T  that  is  a  linear  combination  of  the  single  space  x  =  Xo  operator  a/  given  above.  For 
example,  for  the  case  n  =  2: 
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(o  0  1-2  1, 

Let  A2  be  a  band  diagonal  matrix  with  band  defined  by  a2T,  and 

z  =  A2x. 

Because  the  second-difference  operator  a2T  is  constructed  as  a  linear  combination  of  second-difference 
operators  aj,  the  time  series  z  is  also  stationary  for  the  five  well  known  noise  types  and  linear  frequency 
drift,  and  is  independent  of  both  the  phase  and  normalized  frequency  offsets  for  the  time  series  x.  The 
Allan  variance  for  the  averaging  time  T  =  2to  is  given  by: 

AVAR(t)  =  -^-te[z?] 

3.2  Allan  Variance  for  a  Nonuniformly  Spaced  Time  Series 

Let  (jci,  xs)T  denote  six  consecutive  (uniformly  spaced)  measurements  of  phase  offset  and  consider  the 
situation  that  the  only  values  that  are  available  are  Xu  x4t  and  x6.  Here,  x2,  x3  and  jc5  denote  “placeholders” 
for  missing  measurements.  We  can  derive  a  second-difference  operator  a32T  (the  “3”  and  “2”  represent  the 
3-  and  2-x0  spacing  between  the  measurements  xi  and  x4,  and  x4  and  x6,  respectively)  as  a  linear 
combination  of  single  spaced  T  =  T0  second-difference  operators  a/  that  has  zero  second,  third,  and  fifth 
elements,  as  follows: 
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The  operator  a32  is  determined  uniquely  except  for  a  multiplying  scale  factor  that  is  chosen  so  that  the 
“central”  coefficient  is  -2  (as  for  the  uniformly  spaced  second-difference  operator  a,,1). 

Now  suppose  that  x  is  a  time  series  for  which  the  only  available  measurements  follow  the  above  3-  and  2- 
t0  spacing,  i.e.,  the  available  measurements  are  jq,  x4,  x6,  x9,  xfi,  xi4,  etc.  Let  A32  be  the  matrix 
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z  =  A32x. 
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z  is  calculated  only  in  terms  of  the  available  measurements  in  x.  In  addition,  being  a  linear  combination  of 
single  spaced  x  =  Xo  second-difference  operators  a/,  z  is  stationary  for  the  five  well  known  noise  types  and 
linear  frequency  drift,  and  is  independent  of  both  the  phase  and  normalized  frequency  offsets  for  the  time 
series  x. 

For  this  time  series  we  can  also  construct  a  second-difference  operator  a23T  that  applies  to  measurements 
exhibiting  a  2-  and  3-x0  spacing,  e.g.,  the  measurements  xA,  x6<  and  x9.  Applying  both  a32T  and  a23T 
operators  to  the  available  measurements  in  x  gives  the  more  natural  counterpart  to  the  analysis  for  a 
uniformly  spaced  time  series  where  there  is  “overlap”  between  successive  second-difference  operators. 

An  expression  for  a  general  spacing  of  available  measurements  is  given  as  follows.  Let  (jq,  ...,  x[+p+lIjl 
denote  measurements  of  phase  offset  for  which  the  only  values  that  are  available  are  xu  x[+p,  and  xl+p+l/,  i.e., 
the  spacing  between  available  measurements  is  Ti  =  px o  and  t2  =  qt0.  A  second-difference  operator  apqr, 
that  is  a  linear  combination  of  single  spaced  x  =  x0  second-difference  operators  a/  and  operates  only  on  the 
available  measurements,  is  defined  in  terms  of  these  measurements  by: 


JV-x-2x  -ZZ- 

,  xi  LX\ +P 

p+q  p+q 


Xl+p+q  • 


Now  suppose,  for  a  given  time  series  x,  that  z  contains  all  possible  values  z  computed  by  applying  V  and 
a qp  to  x,  i.e,  from  available  measurements  in  x  separated  by  Ti  =  px0  and  x2  =  <7X0.  Then  we  generalize  the 
definition  of  the  Allan  variance  for  a  uniformly  spaced  time  series  to  be  applicable  to  a  nonuniform  spacing 
as  follows: 

GAVAR^t^-^eIz*  1  r  =  ^£± A 


In  terms  of  its  dependence  on  Xi  and  xt,  this  generalized  Allan  variance  GAVAR{x{,x2)  has  the  following 
properties: 

a)  If  Xi  =  x2  =  x,  then  GA  VAR(xu  x2)  is  identical  to  AVAR(x)  with  an  averaging  time  x. 

b)  The  underlying  time  series  z  on  which  the  GAVAR(xu  x2)  variance  is  based  is  stationary  for  the  five 
well-known  noise  types  and  for  linear  frequency  drift. 

c)  GAVAR(x  1,  x2)  is  independent  of  the  phase  and  normalized  frequency  offsets  for  the  time  series  x. 

Figure  1  compares  values  of  GAVA/?(Xi,  x2)  with  those  of  AVAR(x)  for  WPM,  WFM,  RWFM,  and  linear 
frequency  drift  (LFD).  For  various  choices  of  p  and  q  the  ratio  of  GAVAR(px0,  qx0)  to  AVAR(x),  with 
x  =  x0(p  +  q)/ 2,  is  calculated  and  the  figure  displays  values  of  this  ratio  as  a  function  of  r  =  p/(p  +  q).  As 
expected,  when  p  =  q  (r  =  Vi),  the  ratio  is  unity,  showing  that  the  two  measures  are  identical.  Furthermore, 
for  p  close  to  q,  the  ratio  remains  close  to  unity.  However,  as  r  departs  from  V2  the  ratio  increases  for 
WPM  and  decreases  for  WFM,  RWFM,  and  LFD.  The  curves  for  RWFM  and  LFD  are  very  similar. 
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3.3  Allan  Variance  for  TWSTFT 


A  time  series  of  TWSTFT  measurements  consists  of  measurements  made  on  the  Monday,  Wednesday,  and 
Friday  of  each  week.  For  this  spacing  of  measurements,  we  can  apply  the  second-difference  operators  a 22T 
(=  a2T)  to  the  measurements  on  Monday,  Wednesday,  and  Friday  of  each  week,  a23T  to  the  measurements 
on  Wednesday,  Friday,  and  Monday,  and  a32T  to  the  measurements  on  Friday,  Monday,  and  Wednesday  to 
obtain  a  vector  z.  An  Allan  variance  is  then  calculated  in  terms  of  the  elements  of  z. 

3.4  Modified  Allan  Variance 

In  many  situations,  particularly  in  the  presence  of  WPM,  the  use  of  a  modified  Allan  variance  is  preferred 
to  the  Allan  variance  for  the  characterization  of  TWSTFT  and  other  time  transfer  measurements.  The 
reason  is  that  it  is  not  possible  to  use  the  Allan  variance  to  distinguish  between  the  noise  types  of  WPM 
and  WFM  [7]. 


RATIO  OF  SENSITIVITY  TO  NOISE  TYPES 


♦  WPM 
+  WFM 

♦  RWFM 

♦  LFDM 


Figure  1:  Relative  sensitivity  of  GAVAR  and  AVAR  to  the  five  standard  noise  types  and 
to  linear  frequency  drift. 


The  second-difference  operator  used  as  the  basis  of  the  calculation  of  a  modified  Allan  variance  is 
constructed  as  an  “average”  of  Allan  variance  second-difference  operators.  For  example,  for  a  uniformly 
spaced  time  series  and  averaging  time  x  =  2xo,  we  form 


f1  1  -1 

1 

f  1 

n 

(\ 

0 

-20  1  o' 

U  2 

2 

2) 

u 

2) 

,0 

1 

0  -2  0  1, 

Let  M2  be  a  band  diagonal  matrix  with  band  defined  by  m27,  and 


w  =  M  2\. 
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Then,  the  modified  Allan  variance  MVAR(t)  is  defined  by: 

MVAR(t)  =  ~e[w^\ 

2  t2 

In  a  similar  way,  a  modified  Allan  variance  may  be  defined  for  a  nonuniformly  spaced  time  series. 
Consider  the  case  of  TWSTFT  measurements  made  on  the  Wednesday  (W),  Friday  (F),  and  Monday  (M) 
of  two  consecutive  weeks.  Two  Allan  variance  second-difference  operators  for  the  available  TWSTFT 
measurements  are  as  follows: 
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The  first  of  these  is  for  the  spacing  p  =  5  and  <7  =  4.  The  second  is  for  the  spacing  p  =  q  =  5  and  is  identical 
to  a  second-difference  operator  for  a  uniformly  spaced  time  series  with  n  =  5.  The  second-difference 
operator  used  as  the  basis  of  the  calculation  of  a  modified  Allan  variance  for  the  TWSTFT  measurement  is 
then  given  by  averaging  these  operators,  viz.. 
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Properties  of  a  modified  Allan  variance  for  a  nonuniformly  spaced  time  series  calculated  in  terms  of  such 
second-difference  operators  are  under  investigation. 


4  ESTIMATING  STEPS  IN  TWSTFT  TIME  SERIES 

Delay  steps  of  unknown  magnitude  may  occasionally  be  observed  in  TWSTFT  time  series.  These  are 
usually  due  to  the  replacement  of  a  failed  component  within  the  TWSTFT  instrumentation.  In  some  cases 
it  is  possible  independently  to  determine  the  magnitude  of  a  delay  step,  but  usually,  for  example  when  the 
step  is  the  result  of  changing  a  satellite  transponder,  this  is  not  possible.  We  examine  here  a  method  to 
determine  an  estimate  of  the  delay  step  from  the  time  series  of  TWSTFT  measurements. 

Figure  2  shows  a  time  series  of  (UTC(NPL)  -  UTC(USNO))  TWSTFT  measurements.  Three  delay  steps 
of  unknown  magnitude  are  observed,  at  MJD  51968,  MJD  51984,  and  MJD  52032.  The  first  two  of  these 
steps  are  due  to  unknown  transponder  delay  changes,  while  the  third  is  due  to  a  modem  replacement. 
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4.1  Analysis  Method 


Let  x  =  (xi,  xm)T  be  a  time  series  of  TWSTFT  measurements  corresponding  to  (not  necessarily 
uniformly  spaced)  times  t  =  (6,  ...,  tm)r.  Suppose  delay  steps  in  the  time  series  occur  at  the  three  known 
times  T\,  T2,  and  7>  We  model  these  measurements  in  terms  of  an  nth  order  polynomial  pit,  c)  with 
coefficients  c  =  (ci,  c„)T  and  parameters  d  =  (du  d2,  d2)'  representing  the  (unknown)  delay  steps,  as 
follows: 


UTC(NPL)  -  UTC(USNO)  DELAY  STEPS 


Figure  2:  TWSTFT  measurements  showing  delay  steps  on  MJD  51968,  51984,  and 
52032. 


x,  =  p(ti,c)  +  ei,  t;<T 

X,  =  Pit,  ,c)  +  dx  +  e, ,  7,  <  t,  <  T2 , 

xt  =  p(t, ,c)  +  dl+d2+ei,  T2  < t,  < 73 , 

x,  =  p(f,,c)  +  d,  +d2  +d3  +<?,-,  T^tj. 


(5) 


If  there  were  additional  delay  steps,  the  number  of  parameters  d  and  the  model  (5)  would  be  modified  in  a 
natural  way.  The  measurement  errors  e„  /  =  1,  ...,  m,  are  assumed  to  be  samples  of  WPM  (although  other 
noise  types  and  combinations  of  noise  types  may  be  considered:  see  Section  2).  In  this  case,  estimates  of 
the  parameters  c  and  d  are  obtained  by  solving  the  linear  least-squares  problem: 


m 

minc.d  X*'2- 
1=1 

Equivalently,  the  problem  is  to  find  the  least-squares  solution  b  =  (c1,  d1)1  to  the  linear  system  of  equations 

x  =  Bb  +  e  (6) 

(cf.  equation  (3)),  which  is  overdetermined  provided  m  (the  number  of  measurements)  exceeds  n  +  3  (the 
number  of  polynomial  coefficients  plus  the  number  of  delay  step  parameters).  Here,  the  ith  row  of  the 
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matrix  B  contains  in  its  first  n  columns  the  values  at  t,  of  basis  functions  <J>,(/)  used  in  the  representation  of 
the  polynomial  p(t,  c),  viz., 


p(t,C)  =  Y,Cjtyj(t), 
j=i 


and  in  each  of  its  last  three  columns  the  value  zero  or  one  according  to  whether  the  offset  dj,j  =  1,  2,  3, 
appears  in  the  (corresponding)  ith  model  equation  (5).  One  choice  for  the  basis  functions  <|>y(r)  are  the 
monomials  x!~l,  although  it  is  preferable  to  use  alternative  functions  that  ensure  the  reliability  of  the 
numerical  computations  performed  (see  below). 

Formally,  the  solution  to  this  problem  is  given  by  [2] 

b  =  (flTfl)rVx, 


with  covariance  matrix 


V(b)  =  S2(Brfi)~‘, 

where  s2,  the  root-mean-square  (RMS)  residual  error 
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evaluated  at  the  solution,  estimates  the  variance  of  the  WPM  noise  process.  In  practice,  to  ensure 
reliability  of  the  computed  results,  the  polynomial  p(t,  c)  is  expressed  in  terms  of  Chebyshev  polynomial 
basis  functions  in  a  normalized  variable  [8],  and  the  least-squares  problem  defined  by  (6)  is  solved  using 
matrix  factorization  methods  [9]. 

Using  a  12th  order  (degree  11)  polynomial  to  represent  the  data  shown  in  Figure  2,  the  following  results  are 
obtained: 


d\  =  14.8  ns,  with  standard  uncertainty  1.3  ns, 
d2  =  -28.5  ns,  with  standard  uncertainty  1.6  ns,  and 
di  =  -20.1  ns,  with  standard  uncertainty  1.5  ns. 

In  Figure  3  we  show  the  time  series  of  measurements  together  with  the  fitted  polynomial  curve  following 
the  removal  of  the  unknown  delay  steps. 

Although  polynomials  can  be  effective  for  modeling  time  series  of  TWSTFT  measurements,  they  will  not 
always  be  appropriate.  The  use  of  polynomial  splines  [10]  provides  a  more  general  capability. 
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UTC(NPL)  -  UTC(USNO)  DELAY  STEPS  REMOVED 


Figure  3:  TWSTFT  measurements  corrected  for  unknown  delay  steps. 
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Figure  4:  RMS  residual  error  as  a  function  of  the  polynomial  order. 


4.2  Choosing  the  Order  of  the  Polynomial 

An  important  consideration  for  this  analysis  method  is  the  choice  of  polynomial  order.  Figure  4  shows,  for 
the  example  data  of  Figure  2,  the  RMS  residual  error  s2  as  a  function  of  polynomial  order  n.  The  value 
chosen  for  n  is  the  smallest  value  for  which  the  RMS  residual  error  remains  essentially  constant.  For  the 
example  data  considered  here,  n  is  chosen  to  be  12,  for  which  the  corresponding  RMS  residual  error  is 

1.2  ns. 
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SUMMARY 


In  this  paper  we  have  considered  the  least-squares  analysis  of  measurements  of  the  phase  difference 
between  two  clocks  or  time  scales,  with  the  aim  of  determining  estimates  at  any  epoch  of  the  phase  offset 
between  the  clocks  or  time  scales.  We  have  focused  on  presenting  approaches  developed  to  overcome 
present  operational  constraints  of  TWSTFT,  i.e.,  the  consistently  but  unevenly  spaced  nature  of  the  time 
series  of  measurements,  as  well  as  the  presence  of  occasional  delay  steps  of  unknown  magnitude.  To 
characterize  the  TWSTFT  measurements,  we  have  used  second  difference  statistics  that  are  calculated  only 
in  terms  of  the  available  measurements.  In  order  to  remove  the  unknown  delay  steps,  we  have  modeled  the 
measurements  using  empirical  functions  together  with  parameters  representing  the  delay  steps. 
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